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Summary 


Implicit and explicit spatial differencing techniques with fourth order accuracy have been 
developed. The implicit technique is based on the Pade compact scheme. A Dispersion 
Relation Preserving concept has been incorporated into both of the numerical schemes. 
Two dimensional Euler computation of a spatially-developing free shear flow, with and 
without external excitation, has been performed to demonstrate the capability of numer- 
ical schemes developed. Results are in good agreement with theory and experimental 
observation regarding the growth rate of fluctuating velocity, the convective velocity, and 
the vortex- pairing process. 


Introduction 

Jet noise has emerged as a major issue in the High Speed Civil Transport Program. Noise 
reduction to an acceptable level set by FAR 36 Stage III is a challenging task to achieve. 
In order to accomplish this task, it has become increasingly important to be able to predict 
the sound pressure as a measure of noise level in the aeroacoustic research area. Far field 
noise is generated as a byproduct of the jet flow behind the exhaust nozzle. Therefore, 
understanding of this sound source is crucial. Lighthill [1], through his pioneering acous- 
tic analogy, identified the sound source to be a fluctuating Reynolds stress. Numerous 
other works, both experimental and analytical, have been focused on the flow turbulence 
associated with the sound at a distance. 

Experimental observations show that the sound power emitted from the jet column is 
greatest within 4 or 5 diameters downstream, and then decays rapidly through a transition 
region. This indicates that the initial development of the jet, before it becomes fully 
turbulent, should be clearly understood so that an accurate noise prediction can be made. 
The fully turbulent flow assumption yields analytic results which are far from reality in the 
developing regions. Freymuth [2] observed organized large eddy structures in a separated 
flow of a jet. Brown and Roshko [3] also found large vortical structure in a free shear layer. 
The identity of this large vortical structure is discemable even in the fully turbulent region 
far downstream. Winant and Browand [4] reported that a mechanism of the mixing layer 
growth is an interaction of adjacent large vortices. These investigators have shown that 
the flow in free shear layers such as jet and plane mixing flow is well hehaved and more 
organized than previously thought. Shear flow is dominated by large vortical structures, 
which are very predictable and controllable. These flows, which have been thought to 
be fully turbulent and therefore random and chaotic, have become research subjects with 
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a quite different perspective since the observation of these organized structures. The 
separation of the organized flow entities from the fully random quantities make it possible 
to study turbulent shear flows in a certain deterministic way. 

A turbulent flow field is made up of a range of length scales from the Kolmogorov scale 
to the integral scale. If numerical mesh size can be made fine enough to resolve the 
smallest scales which dissipate the kinetic energy, then direct numerical simulation (DNS) 
is the tool to obtain the entire turbulent flow structure. However, the dissipative scale 
becomes finer as the Reynolds number is increased and the practical hardware limitation 
is rapidly reached. Therefore, the DNS method is limited to simulating only low Reynolds 
number turbulence. For practical computation of higher Reynolds number flows, small 
scale fluctuations can be modeled so that desired large scale eddies can be computed 
directly, while proper dissipation is provided by the small scale eddy model. This approach, 
which is referred to as large eddy simulation (LES), has been successfully employed in many 
flows with practical applications. 

In order to obtain the flow field as the source of sound using DNS or LES the simula- 
tions must be performed using numerical techniques with least distortion and diffusive 
characteristics. The source of numerical diffusion and phase error is known to be mainly 
from the numerical formulation of convective terms. This numerical artifact gets worse 
for high Reynolds number flow simulations. Typically, free shear flows of interest have 
very high Reynolds numbers. Therefore, a higher order accurate numerical scheme which 
meets the previously mentioned requirements is needed. To develop a highly accurate nu- 
merical scheme for this purpose, it will be appropriate to consider only an inviscid flow. 
Fourth order explicit and implicit differencing schemes with a Dispersion Relation Pre- 
serving property are introduced in the context of the numerical formulation of the Euler 
equations. A Pade compact scheme is used in the implicit differencing formulation. A 
plane shear flow generated by two streams of air with different velocities is chosen as an 
example to validate numerical schemes. Subsonic and supersonic results will be shown 
and discussed here together with physical observations of the spatial growth of fluctuating 
velocity, the development of large vortical structure, and the vortex-pairing process. 

Governing Equation 


The variables T, p, x, u, t,p, and e represent dimensionless quantities of temperature, den- 
sity, position, velocity, time, pressure, and total energy per mass. The nondimensionaliza- 

tion is based on reference quantities of T*, p*, l*, u*,p*, and e*. Additional definitions are 
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t* = /*/u*, p* = p*u* , and e* = u* . Then, the non-dimensional equation of state can be 
written as follows : 

pT 


P = 


7 M r 2 


with T = 7(7 — 1 )M? (e — 


u 2 + v 2 


) 


where 7 is the ratio of specific heats, M r = u*/y/jR*T*, and R* is the gas constant. The 
Euler equations in two-dimensional Cartesian coordinates (x, y ) are : 
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s is the source term. The dimensionless speed of sound, c, becomes y/T j M r . Equation (1) 
can be written in a generalized coordinates (£, 77 ) as : 


5Q 5F dG . n t c To 

— ^ 4 \~ —— = S where Q = Jq, S — Js 

dt di dp 

J = x^y v - y^x v = (£ x rjy - ^x) -1 , F = J(^f + ^ y g), G = J(r/ X f + ry y g) 
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Formulation of Difference Scheme 

Basically, two different formulations are presented for a spatial differencing, one is the 
explicit and the other is implicit differencing technique. A Pade compact differencing 
formulation is used for the implicit method developed. Pade compact scheme [5,6] has 
also been previously used in flow computations. All of the methods presented are fourth 
order accurate in space and incorporated with the Dispersion Relation Preserving concept 
proposed by Tam and Webb [7]. 


Explicit Differencing 


The fourth order finite- difference formulation of which is denoted by / , is given by : 



/i+l ~ /t-1 

2 h 


\ _ 1 ( fi± 2 ~ /»- 2 \ 
' 3 ' 4 h ’ 


where the subscript i is the mesh index, h the uniform mesh size. The above expression 
is solely based on the trunction of Taylor series. In unsteady flow computation, the order 
of accuracy is not the only issue to obtain a reliable solution. Unsteady flows by nature 
contain wave trains, which must be well resolved. Numerical difference formulations, con- 
structed solely by truncation of Taylor series, do not guarantee wave preservation even 
when higher order differences are used. The wave becomes severely distorted as the wave 
number, or the frequency in a time space, is increased and eventually a physical inter- 
pretation of the results becomes impossible. Tam and Webb [7] proposed a difference 
scheme, which is called the Dispersion Relation Preserving (hereinafter it will be called 
DRP) scheme, in which the dispersion characteristics follows closely to that of original 
differential form. For example, suppose we need a fourth order explicit difference scheme 
using six neighboring points, then the finite difference formulation of f'{x ) becomes : 
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To be fourth order accurate, the constants b and c are 


and 
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above expression uses four adjacent points, which is solely based on the truncation of 
Taylor series. Equation (3) contains one free parameter, say a, which will be optimized. 
To introduce the DRP scheme, we define the Fourier transform of f(x) and its inverse as : 


K w ) = 

/(*) = 




f(x)e~ iwx dx 
f(w)e iwx dw 


Fourier transforming equation (3) by the above definition gives : 

-v b c ~ 

iK,f(w) = i(a sin k + - sin 2k + - sin 3k)/(u>) 

2 o 

where k is defined to be wh. Notice that with the difference approximation given by (3), 
an input wave number of k is distorted into a different wave number of k which is defined 
to be a sin k + f sin 2k + f sin 3k. It is needed to optimize the free parameter a so that 
the responding wave k depicts the input wave k closely over the wave number space of 
interest. In reference [7] the function to be minimized is chosen as : 


1 = 



( 4 ) 


The limits of integral — ^ and ^ is selected under an assumption that a minimum of four 
mesh intervals are required to resolve a wave. This appears quite plausible because the 
value of I becomes very small over a wide range of a for longer waves, say 8 h, 16 h waves. 
For these longer waves optimization based on DRP concept becomes less meaningful. 


Figure 1 shows the numerical behavior of k against k for different values of a. The straight 
line represents the exact differentiation. It is seen that the curve for a = 1.59853, which 
has been evaluated to be optimum value by reference [7], closely follows the straight line 
for k up to about 0.5, while the sixth order counterpart is valid up to about k of 0.4. A 
conventional fourth order differencing, which is obtained at a = |, restricts its validity 
only to a long wave range, say k < 0.3. There is no doubt that the second order accurate 
scheme is severely limited to resolving only long waves. 


Pade Compact Differencing Scheme 


The Pade compact difference scheme used here is an implicit method in space. This scheme 
uses a symmetric numerical stencil in a compact form so that only a tridiagonal or penta- 
diagonal matrix inversion is needed. Even for higher order formulations the compactness 
is still maintained. Following general form of the Pade difference scheme can be used to 
obtain a difference form of /'. 


• • • + + f'i + + Ol2fi + 2 + 
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In the present formulation we use the following form to assure that only the tridiagonal 
matrix inversion is required to obtain a solution for 
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With a = 


and b = ^ 


t this formulation is forth order accurate, and a ia a free 

parameter. Equation (5) can be reduced to sixth order accuracy if a is set to |. It is 
straightforward to make the Pade formulation (5) DRP based since it has a built-in free 
parameter a to be optimized. Using the same definition of Fourier transform as used in 
the explicit differencing case, we can formulate the wave number relation as : 


K = 


a sin k 4- |sin2/c 
1 + 2a cos k 


The optimum value of a is computed to be 0.35619 when the integral I, defined in (4), is 
minimized. Figure 2 illustrates the wave number relation between k the input wave and 
k the response wave. As was the case with the explicit difference scheme, better wave 
behavior is observed for the fourth order scheme with the optimum a than for the sixth 
order differencing with a = |. The curve for a = j is of the most compact fourth order 
scheme. Figure 3 is the comparison made among the six cases. The most compact fourth 
order scheme has broader range of wave number than the optimum case of the explicit 
scheme. However, for a computation which employs a large number of mesh points such 
as in a three dimensional flow simulation, the optimized explicit scheme can be a good 
alternative to the implicit differencing if only a moderate number of mesh points are 
added more for the resolution of all waves under consideration. 


Comments on Finite Volume Scheme 


Finite volume formulations have been favorably used in the CFD community because it 
mimics an integral form of physical equation on a computation cell. The surface values, 
known as fluxes, are of primary concern. The finite volume formulation of the flow equation 
is conservative since the fluxes are cancelled out when a global summation is carried out 
over all the computation cells. For example, the one dimensional equation, -57 + = 0, 

can be integrated over a space interval of h to give : 



udx + F i+ 1 




1 

2 


= 0 


where F is a function of u, h uniform mesh size, and the subscripts i - \ and 1 + \ denote 
the left and right cell boundaries. A second order spatial accuracy can be achieved by 
representing the integral, J u dx, by Uih and the fluxes F 1+ i and F } _i to be \{F i+l + Fi) 
and |(F,_i + Fi). This leads to the symmetric finite difference formulation of : 

— + — - + Fi+1 — i? *~ 1 + 0 ( h 2 ) 

dt dxdt 2 h K ' 
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For higher order spatial accuracy, the finite volume formulation is proposed as the following 

[ + * ,dui dF » lV'' 9 . ^ 

L { -si + ^ ix = h E 

2 k= — N 

Expanding F t _ i and F i+ i about the point x, in the Taylor series and then using F'" = 
—du"/dt, F'"" = —du""/dt (prime denotes a derivative with respect to the spatial variable 
x) relationships, we can obtain the following for fourth order accuracy. 


du dF d , 1 22 1 

m + ^ ~ S ( 24“-‘ + 2i Ui + 24“ i+l) + 


F *i - F <- 


+ 0(/i 4 ) 
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The quantity in the parenthesis of the time term is a cell average value of u and the fluxes 
are evaluated in a symmetric manner as: 



9 (Fj+i + Fj) — (Fj+ 2 + Fj-i) 
16 


and F.i 
1 2 


9(F i . 1 +F i )-(F i . 2 + F i+1 ) 
16 


Equation (6) is a time implicit form which calls for a matrix inversion. A multi-dimensional 
extension of the finite volume formulation becomes increasingly complicated. 


As a coding practice in the present work, a finite volume like formulation is used for the 
spatial discretization through a definition as : 

, fi+% ~ h-\ 

J h 



The value / can be treated as a cell interface value only by analogy, but it is not the 
physical value at the designed accuracy, (up to the second order accuracy, / represents a 
physical quantity at the interface) For the explicit differencing, we propose the following : 


r a fi + /•- 1 , fi+1 + fi-2 . . fi+2 + fi-z 

Ji- 1 = Ai b A 2 b A z 


Using the definition (7) for equation (3), we arrive at the following constants : 

, be be. c 

Al=a+ 2 + 3’ 2 = 2 + 3’ Aa = 3 

For the Pade compact scheme, we also propose the following formulation, which is similar 
to equation (5), to obtain an interface value at the point i — \- 


£ i 2 i £ 7j /* f" fi 1 i d f i+ 1 f" / i — 2 

+ Ji-k + a fi+h = B 1 o + B * ^ 


( 8 ) 


Rewriting a similar equation centered at i + h and then subtracting the original equation 
(8) from it gives : 


a 


(fi-h -fi-O + fi+i ~ 1 ~fi+0 = B 2 


fi+2 ~ fi- 




fi - 1 
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( 9 ) 


By using the definition of equation (7), the above equation reduces to 
+ /,' + «/!» = + (B, - B 2 ) 

By inspecting equation (9) against equation (5), the constants are found to be B\ — 8a 6 +7 , 
and Bo — 4o 6 ~ J . For o = i only two neighboring two points are used, which has the most 
compact form in the right hand side of equation (9). 

Using F and G defined in a similar method as the /, we can spatially discretize the Euler 
equation (2) by setting A£ = At? = 1 as : 


5Q 

dt 


+ F i+ij -F.-ij +G ij+ i 



S 


The indices i and j are for the £ and ?? directions, respectively. 


Time Advancing Scheme 


The four-stage Runge-Kutta technique [8] is adopted for the explicit formulation in time. 
To obtain new flow variables at t — (n + l)At from known data at t = nAt equation (2) 
is used to advance the solution in time as follows : 

- Q n = ati AfW (0) 
q( 2 ) - Q n = a 2 AfW (1) 
qU) _ Q n = a 3 AtW (2) 

Q( 4 ) - Q n = a 4 AtW (3) +D 

where W denotes S — 3F /3£ — dG/di] evaluated at the fc-th stage. The stage 0 and 
4 are at the time, nAt and (n + l)Af . The parameters, Oi, 02 , 03 , 0:4 , are given to be 
- - 1. This time difference is second order accurate and the intermediate variables are 

not stored at every stage. 

A numerical dissipation term D is added during the fourth stage to enhance the numerical 
stability. The dissipation term is introduced to be of sixth order so that our fourth order 
accuracy remains intact. 

n- t (— q I d 6q l 

D We d£ 6 + dr? 6 ' 

where u e is a constant and d 6 q/9£ 6 is given by : 

= 15(q,-+ij + q^-ij) - 6(q !+2j + q t - 2 >) + (q>+3j + q,-3j) - 20q t y 

The derivative <9 6 q j drf in the r? direction is obtained in the similar manner. This numerical 
dissipation is applied to all internal points. 
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Boundary Condition 


The boundary treatment considered here is a combination of characteristic and algebraic 
boundary conditions. The characteristic boundary condition solves the governing equation 
in a characteristic form in each coodinate direction and the algebraic boundary condition 
is the given boundary conditions such as temperature, total temperature, velocity, etc. 
Equation (2) is written in a non-conservative form as: 

da dq „ dq , . dF 5G 

dt di dr] dq dq 

For the f direction, equation (2) with a transformation dq = R^dq becomes 

| +R( -. AR( | + r-. B r s | = o 

If we construct the matrix such that the eigenvectors of the matrix A constitute 
its col umns then the matrix by a similarity transform becomes a diagonal matrix, whose 
entries are the eigenvalues of A such that 

R^AR^ = A^ = diag(U, 17, U 4- a^, U — a 

where = c yjil + • Here, c is the speed of sound defined to be y/T /M r , and U = 

£_ x u + £ y v. In the same way, for the 77 direction the diagonal matrix becomes : 

R~ 1 BR 1} =* A, = diag(V, V, V + a„ V - a v ) 

where a v = c yjr\\ + rjy , and V = rj x u+r]yV. The characteristic equations are then rewritten 
in each coordinate direction to be: 


t -idq 
« dt 


R* = 


dq 

di 

♦? 


and R-‘| + A,R-| + 
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v4c \Fic 
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PK y 
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—pk x 
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J d^ 
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V^rr p(uk y -vk x ) 




V--Z-8 
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'-ho 

r—H 
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p 

P 
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-e , $ 

•\/2p V5 pc 

, (I- 7 )* 

V5p \/2pc 

1 (i-7)w 
y/2 p \/2pc 

7-1 

V2pc 


V \/2p V^pc 

(l-7)u 
\/2p V2pc 

Z*JL (1-Tf).» 
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\/2 pc 
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9 = K X U + KyU, 


7-1 


(u 2 +t> 2 ), 


V^( 7 -l)c n/2( 7 -1) 


Hedstrom [9] proposed that the non-reflecting boundary condition can be constructed by 
setting any eigenvalue, which is the element of A, to be zero, if a wave is incoming towards 
the computational domain. Thompson [10] further extends this method and Lele [5] used 
the same boundary condition in the simulation of a free shear flow. 


Application 

In the present study Euler computations of plane shear flow are presented to validate 
the numerical method developed. For this validation, computional results of spatially 
developing flow are examined against linear theories, experimental data, and flow visual- 
izations. The initial shear is given as the axial velocity distribution at x = 0, where x is 
the downstream variable. Initial vorticity thickness of the shear SijJq, which is defined as 
( u* - u*)/(du*/dy*) max , is chosen as the reference length l*. Quantities with asterisks are 
dimensional, subscripts / and s denote faster and slower speeds, respectively. The maxi- 
mum derivative is taken at the inflection point at y = 0. Then, the initial shear profile is 
riven as : 

6 Au* 

u = \ 3 tanh(2y) + 1, A 3 = (10) 

where Au* = u* f -u* s , u* = (u* f +u*)/ 2, and the reference velocity in this case is the average 
speed of two streams. The momentum thickness, 9 , which is used as a reference length scale 
in both theoretical and experimental analysis, becomes o/4 for the shear given by (10). 
All the flow computations is carried out by assuming that mixing of two streams occurs at 
an uniform temperature. The other flow parameters chosen are: T* = 298°K, 1*= 0.00254 
m, and p* = 1.1839 Kg/m 3 . 

Figure 4 illustrates two grids used. Grid I (240 x 120) and Grid II (300 x 160) are stretched 
in both y and x directions. In Grid I mesh size Ax varies between 0.4 and 1.568 with 
0.15 < A y < 0.4 at the inlet and 0.333 < A y < 0.864 at the exit plane. The inlet extends 
between -12.43 < y < 12.43 and the exit plane between -32.7 < y < 32.7, and x reaches 
196.5. In Grid II mesh size Ax varies between 0.516 and 0.84. Inlet and exit have the same 
height of -25.8 < y < +25.8 with Ay min = 0.12 and A y max = 0.576. The downstream 
extends to x = 206. 


Subsonic flows are simulated with and without external excitation on the initial shear. 
Mach numbers of the two streams considered here are 0.6 and 0.3, in which u*=155.768 
m/s, M r = 0.45, and p* is 2.873 K Pa. At the inlet, u defined by (10), u = 0, T = 1 
are given and one characteristic equation is solved for the outgoing acoustic wave. At the 
exit, the characteristic equations are solved for the three outgoing waves. For the other 
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incoming acoustic wave, the characteristic equation is also imposed with its incoming wave 
velocity set to zero. On the side boundaries, all the characteristic equations axe computed 
by setting the wave velocity to zero if the wave is incoming towards the computational 
domain. The flow case with no excitation, which is naturally acquired flow, is simulated 
first. A forcing is then given later on the v velocity at the inlet from the natural flow 
solution to excite the flow. 

Figures 5-(a,b,c) show the contour plots of Mach number, vorticity, and the static pressure 
at an instant in time. A large vortical structure is clearly seen in the Mach number and the 
vorticity contour plots as shown by Figures 5-(a) and (b). The static pressure contour plot 
of Figure 5-(c) supports existence of this vortical structure by aligning its local mimima 
with the centers of the vortices depicted in the Mach number and the vorticity contour 
plots. This train of local pressure minima and maxima is the basis of analytic estimation of 
the convective velocity at which the large vortical structure travels. The vorticity contour 
plot is chosen hereafter as an appropriate visual device to examine the large flow structure. 
The vorticity contour plots in Figures 6 through 9 show a spatial evolution of shear flow 
with time. The convective velocity can be found from the time history of the vorticity 
contours. The simulation using the Pade differencing with a = 0.3 and 0.35619 predicted 
the convective velocity, u c , of 0.99. In the explicit computation u c has a value between 0.96 
and 1.02. From the four cases shown in Figures 6 through 9, it can be concluded that the 
convective velocity is about unity, which is the average velocity of the two streams. Brown 
and Roshko [3] found that the large eddies travel at a constant velocity, which is the average 
speed of the two streams, and this velocity is independent of size or location of the eddy. 
Papamoschou and Roshko [11] derived the same result. This convective velocity is identical 
with the computed value and the present results also confirm an earlier observation that 
the convective velocity is independent of eddy location and size. The only exception is 
that the convective velocity varies when two successive vortices merge into one, which is 
often referred to as a vortex pairing. In this pairing process, the vortex which travels 
behind picks up speed to catch up with another vortex proceeding ahead, then the two 
slide on each other and coalesce. After the completion of vortex pairing the convective 
velocity becomes constant again until another pairing occurs further downstream. This 
vortex pairing process is illustrated in the vorticity contour plots. Figure 10 shows the 
trajectories of vortices, which shows the birth, pairing, convecting, and further pairing of 
vortices. Parallel slopes on the x — t diagram, which is the convective velocity, are clearly 
seen to be constant regardless of the vortex location. 

An unforced shear flow exhibits a distributed spectrum in its fluctuating quantity. Figure 
11 shows this for the u velocity spectra at y = 0 and various downstream locations. A 
peak in the spectrum is clearly seen in this unforced flow. The frequency for the peak in 
the spectrum is referred to as the most preferred frequency. The Strouhal number, St. 
is used as a dimensionless frequency and defined to be f*0/u* with f* the frequency in 
Hz, 6 the initial momentum thickness. In the early flow development, the u spectrum 
has its peak at St = 0.021. This persists to some location downstream, then eventually, 
the lower frequency of St = 0.0086 becomes dominant far downstream. This agrees with 
the general observation that lower frequency wave, whose wave length is longer, survives 
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further downstream than a higher frequency wave. Linear theories by Michalke [12], and 
Monkewitz and Huerre [13] show that the most amplified frequency is about St = 0.033. 
In jet flow, Freymuth [2] found the growth rate observed experimentally agrees very well 
with the theory for St up to about 0.024, and remains constant for St > 0.024. This value 
is very close to the computed St of 0.021. In light of those results, it appears that the most 
preferred frequency, which is acquired naturally, has a lower value than the most amplified 
frequency predicted by linear theory. 

It is known from linear theory that a free shear flow is unstable under presence of a dis- 
turbance. The disturbance introduced grows exponentially either in space or in time. The 
initial shear layer undergoes wave propagation in the presence of a background disturbance, 
which is due to the truncation error caused by a discretization of the flow equations. Trav- 
eling instability waves get intensified to a certain level of magnitude. In Figures 6-9, the 
shear layer appears wavy, yet it is well connected up to x of 60 to 80 depending upon 
time. The shear layer, then, rolls up into a discrete vortex and convects downstream. In 
the early flow development stage before the vortex roll-up, there is a region where linear 
theory is valid. This linear region can be verified by checking for exponential growth in the 
magnitude of fluctuating u velocity. The root-meam-square of fluctuation of u, denoted 
by uj. ms , is used as a measure of the magnitude. Figures 12-(a),(b), and (c) show the 
region in which a logarithmic scale of grows linearly for various y locations. Sharp 

increase near x — 0 should be neglected since the boundary condition imposed at the inlet 
influences the flow solution in this region. It appears that the slope vanes also with y and 
the maximum slope is shown on the slower fluid side. A growth rate s is defined to be 
dilnu'-nn^/dx. Then, values of s for the unforced shear flow shown in Figure 12-(a) are 
0.096, 0.102, 0.106, 0.105, 0.1, 0.099, and 0.096 at y=-4, -3, -2, -1, -0.06, 1, 2, respectively. 
Figure 12-(b) presents curves for when forcing is applied at its most preferred upstream fre- 
quency at St of 0.021. Values of the slope are the same as those for unforced flow, because 
natural shear flow is already dominated by the wave of that frequency. However, the linear 
region near the inflection point becomes hardly discemable. This can be an indication that 
nonlinearity takes place earlier than the unforced case. With lower frequency forcing, s 
has smaller values of 0.087, 0.096, 0.098 at y=-3, -2, -1, respectively. It is difficult to find 
the constant value of s near the inflecton point and in the faster fluid territory. It can be 
interpreted that the two waves, one of forcing frequency and the other of upstream most 
preferred frequency, interact nonlinearly. However, it appears in common that grows 

up to a certain magnitude, and then levels off. The uj^s varies in a mildly oscillatory 
manner further downstream, which is not shown in the figure. In Figure 12-(a) and (c), 
u rms becomes magnified up to about 10 %, and 13 % in Figure 12-(b). 

The uj-ms discussed above is not for a single frequency, but for all the frequency content. To 
observe the growth behavior of individual waves a Fourier decomposition is needed. This 
has been done for the u velocity at y = 0 at various downstream distances. The growth 
rates s computed are 0.065, 0.108, 0.12 for St=0.0086, 0.021, 0.028, respectively. The 
exponential growth of the low frequency wave at 5^=0.0086 is observed about 65 < x < 75, 
which is sufficient downstream where the vortex roll-up can occur. Linear wave behavior 
is found in the region of x < 36 and x < 25 for St of 0.021 and 0.028, respectively. Low 
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frequency expansion given by Monkewitz and Huerre [13] yields the growth rate s(St) at 
X a of | to be 0.578, 0.111, 0.122, in the same order as above. 

The profile of tt^g is illustrated in Figure 13. In the early development region kinks 
appear on both sides of the peak near the center. The kink on the slower side of fluid 
is very sharp. This phenomenon has been reported in both the experimental literature 
[2,4,15] and the theoretical literature [12]. Figure 13-(a) is for the unforced flow. The 
u rms profile first develops into double humps, then evolves into a single peak further 
downstream. Figure 13-(b) also shows the same type of kink in the initial development, 
but quite a different shape forms downstream. With excitation at the preferred frequency 
a deep valley appears in the center area. This is similar to the phenomenon Browand 
[15] found when the shear flow is excited at near the most preferred frequency. As flow 
proceeds further downstream the profile evolves to a single peak as in the previous case. 
As shown in Figure 13-(c), «rms under a lower frequency forcing, which is away from the 
most preferred upstream frequency, takes on a character similar to the unforced flow. No 
deep valley is seen in the middle of the flow evolution. 

Forcing is introduced to the initial shear to see the more organized flow structure. The 
transverse velocity v is given an oscillation of cg(y) sin u>t , where g is a Gaussian distribu- 
tion. For a subsonic mixing of two streams of M=0.6 and 0.3 e is given to be 0.01, and the 
forcing is given at the most preferred frequency of St =0.021. The angular frequency u? is 
which is 5.2 K Hz for the flow condition given. Figure 14 is the vorticity contours 
at consecutive time periods 192 n seconds apart. As shown in Figures 6-9, vortices do not 
appear to be periodic in the unforced flow, the spacing between two consecutive vortices 
is irregular, and vortex roll-up and merging of two or even three vortices occurs randomly. 
However, in the presence of forcing, the initial vortex pairing process is suppressed, with 
the process resumed at some downstream location. Vortices are orderly, keeping the same 
distance apart. This is similar to the situation that Ho and Huang [14] observed. As 
shown in Figure 14, forcing generates a very organized flow structure. This figure also 
shows that the first five vortices appear frozen, in other words the flow up to about five 
vortices downstream is stationary. That is the precise counterpart of the flow visualization 
made in the laboratory using conditional sampling technique, which can be obtained by 
synchronizing the strobe speed with the forcing frequency. Vortex pairing occurs down- 
stream, consequently, the flow appears to be unsteady in shedding vortices. This vorticity 
contour also gives the convective velocity to be 0.99, which is the same value found for the 
unexcited case. 

A forced supersonic free shear flow is also computed with M=1.6 and 1.2. All of the velocity, 
temperature, and density are given at the inlet. All of the characteristic equations are 
imposed at the exit as the boundary conditions. Unlike subsonic flow, it is very difficult (at 
least in the present computation) to simulate the supersonic shear flow without excitation. 
This indicates that the nature of free shear flow generated by two supersonic streams is 
far more stable than the subsonic flow. External excitation of e = 0.05 is imposed on the v 
velocity at the inlet in a similar manner as used in subsonic flow. The forcing frequency u> is 
taken arbitrarily to be Figure 15 shows the flow structure at every two time periods, 
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which is at 128.9 n seconds apart. The supersonic computation has been performed on 
Grid I. It is interesting to observe that no vortex pairing occurs. Many experiments have 
revealed that the free shear of two supersonic streams exhibits poor mixing characteristics, 
which is an indication of a very stable flow situation. All of the vorticity contour appears 
identical because conditional sampling is made on the flow field which is free of vortex 
paring processes. Again, the convective velocity computed from Figure 15 is 0.99. This 
value is in agreement with Papamoschou and Roshko [11]. From these observation of the 
excited shear flows it can be concluded that the forcing has no influence on the convective 
velocity. 


Future Direction Towards Aeroacoustics 


The fourth order numerical schemes developed in this work present flow solutions which 
exhibit important flow phenomena and agrees very well with published experiment and 
theory. As mentioned briefly in the introduction the flow field so simulated will be used 
as a source of the sound. There are two main ways to predict far-field noise from the 
n um erical flow solution. The first approach is the acoustic analogy by Lighthill [1]. In this 
method, a numerical solution of the wave equation for a far-field is not required since the 
integral form of the exact solution is used, which is called the retarded potential solution. 
The computed flow solution is fed into the integrand and is to be integrated over the 
entire flow domain. Since the integrand is evaluated at different retarded times, a time 
interpolation is needed. Therefore, flow solutions must stored at several time steps over 
the entire flow regime. 

A second way is to solve the wave equation directly, which can be done parallel to the flow 
field computation. A flow boundary is imposed sufficiently far from the center of the free 
shear layer, so that acoustic assumption is deemed plausible there. The wave equation 
is solved from a surface near the flow boundary with the pressure or density distribution 
given on it. The mesh in the wave equation computation is much larger than that used 
in the flow computation and becomes stretched greatly. The schematic is drawn in Figure 
16. An accurate wave equation solver is an essential tool for this purpose. A higher order 
DRP scheme seems to be a good choice to solve the wave equation. 

It is acknowledged that Dr S.T. Yu's work was a benefit to the present development of 
numerical methods. 
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Figure 6.— Vorticity contour at every 200 At computed by Pade scheme with optimum a. 
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Figure 7. — Vorticity contour at every 200 Af computed by Pade scheme with a = 0.3. 
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Figure 8.—Vorticity contour at every 200 A t computed by DRP explicit scheme with a = 1 .59853. 
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Figure 9.— Vortlclty contour at every 200 At computed by explicit scheme with a = 4/3. 
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Figure 1 4.— Vortlclty contours of forced subsonic flow at every period. 
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Figure 1 5.— Vorticity contours of forced supersonic flow at every two period. 
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Figure 1 6. — Schematics of flow and wave computation 
regions. 
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